#include "cpgplot.h"
#include "fitsio.h"

#include <iostream>
#include <fstream>
#include <unistd.h>
#include <vector>
#include <cstring>
#include <algorithm>

#include <stdio.h>
#include <stdlib.h>
#include <cmath>
#include <time.h>

#include "posthack.h"

#ifndef NO_LUA
#include "luaclass.h"
#endif

#include "eclipse_types.h"

#ifndef EXIT_FAILURE
#define EXIT_FAILURE 1
#endif
#ifndef EXIT_SUCCESS
#define EXIT_SUCCESS 0
#endif

extern float ContourSmoothing;
extern float ContourStart;
extern float ContourSampling;

enum palettetype {GRAYSCALE, RAINBOW, HEAT, IRAF, AIPS, RED, GREEN, BLUE, BB, BB2, RAINBOW2, RAINBOW3};

enum palettescaling {LINEAR, SQRT, SQUARE, CBRT, CUBE, LOG, ALOG, HISTEQ, POW, ASINH, SINH, TANH, ATANH, LUPTON};
float powpar = 1.0;

void Draw ();
void palett (palettetype, palettescaling, float, float);
void setvp ();
void LoadFile (const char *, float **);
void LoadFitsFile (fitsfile *, float **);
void LoadContours (const char *);
void Resample (float, float **, bool);
bool ParseColor (const char *s, float *r, float *g, float *b);
void TriangularSmooth (float **);
void FileAll (const char* name, float **buf); 
void FluxPerc (float *, float*);

void Help (void);
void DrawWedge (const char *side, float disp, float width, 
		float fg, float bg, const char *label);

std::vector<float> contour_levels;

std::vector<int> cont_colors;
std::vector<int> cont_styles;
std::vector<float> cont_values;

int linewidth;
float minperc = 0.05;
float maxperc = 0.9975;

int nRGB = 0;
float * RGB[3] = {NULL, NULL, NULL};
float * RGB_val = NULL;
float gamm = 1.0;

int sizex=64,sizey=64;
int oldsizex=64, oldsizey=64;
int smaxx = -65536, smaxy = -65536;
int sminx = -65536, sminy = -65536;
float bright = 0.5;
float contra = 1.0;
float smax = -1e10;
float smin = -1e10;
float *data = NULL;
float *cont;
int xmin=65536, xmax=-65536,ymin=65536, ymax=-65536;

float gflmin=-1e10;
float gflmax=-1e10;
float bflmin=-1e10;
float bflmax=-1e10;

int window_size = 0;

bool wedge = true;
bool reverse = false;
bool invert = false;
bool numbers = true;
bool whitebcg = false;
bool setbcgcolor = false;
bool trismooth = false;
bool noxnumbers = false;
bool noynumbers = false;
bool flipx = false;

float bcg_r = 1.0;
float bcg_g = 1.0;
float bcg_b = 1.0;
bool wedgeright=false;

bool verbose = false;
const char *title = NULL;
const char *output = NULL;
const char *xlabel = NULL;
const char *ylabel = NULL;
bool contour = false;
const char *wedgelabel = "";

float flmin=1e20, flmax=-1e20;

float pixelscale_x = 1.0;
float pixelscale_y = 1.0;
bool pixelscale = false;
float xref1 = 0.0;
float yref1 = 0.0;
float xref = 0.0;
float yref = 0.0;
float xrefmin = 0.0;
float yrafmin = 0.0;
float sigma=-1.0;
int region=1;
int iter=1;
float zoom = -1.0;
float c_lupton = 5.0;

float ticklength = 1.0, labelsize = 1.0;
int tickthick = 3, labelthick = 2;

float lgcon = 0.01;

bool axis = true;
bool ticks = true;
bool perc = false;

palettetype ptype = RAINBOW;
palettescaling pscale = LINEAR;

const char *gimage = NULL;
const char *bimage = NULL;

#ifndef NO_LUA
const char *luafile= NULL;
const char *luastring = NULL;
Lua lua;
#endif


int main( int argc, char ** argv )
{
  int n;
  bool postscript = false;
  bool postscripthack = true;

  if (argc == 1) {
    Help ();
    return (EXIT_SUCCESS);
  }

  char buff[64];
  time_t time_t_timestamp; 
  time(&time_t_timestamp);
  char *s_timestamp = ctime (&time_t_timestamp);
  int i=0;
  while (s_timestamp[i] != 0 && s_timestamp[i] != '\n') {
    buff[i] = s_timestamp[i];
    i++;
  }
  buff[i++] = ':'; buff[i]=0;

  FILE *log = fopen (".colormap", "a");
  if (log != NULL) {
    fprintf (log, "%s ", buff);

    for (int i=0; i<argc; i++)
      fprintf (log, "%s ", argv[i]);
    fprintf (log, "\n");
    fclose (log);
  }
  else
    printf ("ERROR: %m\n");

  while ( (n=getopt (argc, argv, "v:f:z:cg:hm:p:a:b:t:o:s:x:y:l:")) != -1) {
    char * temp;
    float tempvalf;
    int tempvali;
    //printf ("%c %s\n", n, optarg);
    switch (n) {
    case 'f':
      if (strcmp(optarg, "noaxis") == 0)
	axis = false;
      else if (strcmp(optarg, "nonumbers") == 0)
	numbers = false;
      else if (strcmp(optarg, "verbose") == 0)
	verbose = true;
      else if (strcmp(optarg, "wedgeright") == 0)
	wedgeright = true;
      else if (strcmp(optarg, "noticks") == 0)
	ticks = false;
      else if (strcmp(optarg, "reverse") == 0)
	reverse = true;
      else if (strcmp(optarg, "nowedge") == 0)
	wedge = false;
      else if (strcmp(optarg, "rticks") == 0)
	invert = true;
      else if (strcmp(optarg, "noxnumbers") == 0)
	noxnumbers = true;
      else if (strcmp(optarg, "noynumbers") == 0)
	noynumbers = true;
      else if (strcmp(optarg, "flipx") == 0)
        flipx = true;
      else if (strcmp(optarg, "trismooth")==0)
	trismooth = true;
      else if (strcmp(optarg, "whitebcg") == 0) {
	printf ("WARNING: -f whitebcg is deprecated!\n");
	printf ("use -v bcgcolor=#FFFFFF instead\n");
	setbcgcolor = ParseColor ("#FFFFFF", &bcg_r, &bcg_g, &bcg_b);
      }
      else if (strcmp(optarg, "24bit")==0) {
	printf ("This option is not needed, as 24bit colours are now used as default\n");
	printf ("Use \"-f 8bit\" if you wish to use old system\n");
	postscripthack = true;
      }
      else if (strcmp(optarg, "8bit")==0) 
	postscripthack = false;
      else if (strcmp(optarg, "perc")==0) 
	perc = true;
      else
	std::cout << "Unknown option for -f: " << optarg << std::endl;
      break;
    case 'v':
      temp = strstr (optarg, "=");
      if (temp == NULL)
	std::cout << "Invalid variable format" << std::endl;
      else if (strlen(temp) <= 1)
	std::cout << "Invalid variable format" << std::endl;
      else {
	tempvalf = atof (&temp[1]);
	tempvali = atoi (&temp[1]);
	if (strncmp(optarg, "contrast=", strlen("contrast="))==0)
	  contra = tempvalf;
	else if (strncmp(optarg, "brightness=", strlen("brightness="))==0)
	  bright = tempvalf;
	else if (strncmp(optarg, "gamma=", strlen("gamma="))==0)
	  gamm = tempvalf;
	else if (strncmp(optarg, "contoursigma=", strlen("contoursigma="))==0)
	  ContourSmoothing = tempvalf;
	else if (strncmp(optarg, "contourstart=", strlen("contourstart="))==0)
	  ContourStart = tempvalf;
	else if (strncmp(optarg, "contoursampling=", strlen("contoursampling="))==0)
	  ContourSampling = tempvalf;
	else if (strncmp(optarg, "tick_length=", strlen("tick_length="))==0)
	  ticklength = tempvalf;
	else if (strncmp(optarg, "tick_thickness=", strlen("tick_thickness="))==0)
	  tickthick = tempvali;
	else if (strncmp(optarg, "label_size=", strlen("label_size="))==0)
	  labelsize = tempvalf;
	else if (strncmp(optarg, "label_thickness=", 
			 strlen("label_thickness="))==0)
	  labelthick = tempvali;
	else if (strncmp(optarg, "wedge_label=", 
			 strlen("wedge_label="))==0)
	  wedgelabel = &temp[1];
	else if (strncmp(optarg, "x2=", strlen("x2="))==0)
	  smaxx = tempvali;
	else if (strncmp(optarg, "x1=", strlen("x1="))==0)
	  sminx = tempvali;
	else if (strncmp(optarg, "y2=", strlen("y2="))==0)
	  smaxy = tempvali;
	else if (strncmp(optarg, "y1=", strlen("y1="))==0)
	  sminy = tempvali;
	else if (strncmp(optarg, "cx=", strlen("cx="))==0)
	  xref1 = tempvalf;
	else if (strncmp(optarg, "cy=", strlen("cy="))==0)
	  yref1 = tempvalf;
	else if (strncmp(optarg, "region=", strlen("region="))==0)
	  region = tempvali;
	else if (strncmp(optarg, "sigma=", strlen("sigma="))==0) 
	  sigma = tempvalf;
	else if (strncmp(optarg, "iter=", strlen("iter="))==0) 
	  iter = tempvali;
#ifndef NO_LUA
	else if (strncmp(optarg, "luafile=", strlen("luafile="))==0) 
	  luafile = &temp[1];
	else if (strncmp(optarg, "luastring=", strlen("luastring="))==0) 
	  luastring = &temp[1];
#endif
	else if (strncmp(optarg, "scale=", strlen("scale="))==0) 
	  zoom = tempvalf;
	else if (strncmp(optarg, "gimage=", strlen("gimage="))==0) 
	  gimage = &temp[1];
	else if (strncmp(optarg, "bimage=", strlen("bimage="))==0) 
	  bimage = &temp[1];
	else if (strncmp(optarg, "gflmin=", strlen("gflmin="))==0) 
	  gflmin = tempvalf;
	else if (strncmp(optarg, "bflmin=", strlen("bflmin="))==0) 
	  bflmin = tempvalf;
	else if (strncmp(optarg, "gflmax=", strlen("gflmax="))==0) 
	  gflmax = tempvalf;
	else if (strncmp(optarg, "bflmax=", strlen("bflmax="))==0) 
	  bflmax = tempvalf;
	else if (strncmp(optarg, "bcgcolor=", strlen("bcgcolor="))==0) 
	  setbcgcolor = ParseColor (&temp[1], &bcg_r, &bcg_g, &bcg_b);
	else
	  std::cout << "Unidentified variable" << std::endl;
      }
      break;
    case 'x':
      xref1 = atof (optarg);
      break;
    case 'y':
      yref1 = atof (optarg);
      break;
    case 'c':
      printf ("WARNING: using contour files is deprecated!\n");
      printf ("         use lua file and Contour ({list}); instead!\n");
      printf ("         the contourfile will NOT be loaded, if it's specified\n");
      contour = true;
      break;
    case 'g':
      if (strcmp(optarg, "linear") == 0) 
	pscale = LINEAR;
      else if (strncmp(optarg, "log", strlen("log")) == 0) {
	pscale = LOG;
	if (strlen(optarg)>3)
	  lgcon = atof (&optarg[3]);
      }
      else if (strcmp(optarg, "alog") == 0)
	pscale = ALOG;
      else if (strcmp(optarg, "sqrt") == 0)
	pscale = SQRT;
      else if (strcmp(optarg, "square") == 0)
	pscale = SQUARE;
      else if (strcmp(optarg, "cbrt") == 0)
	pscale = CBRT;
      else if (strcmp(optarg, "cube") == 0)
	pscale = CUBE;
      else if (strcmp(optarg, "asinh") == 0)
	pscale = ASINH;
      else if (strcmp(optarg, "sinh") == 0)
	pscale = SINH;
      else if (strcmp(optarg, "tanh") == 0)
	pscale = TANH;
      else if (strcmp(optarg, "atanh") == 0)
	pscale = ATANH;
      else if (strcmp(optarg, "histeq") == 0)
	pscale = HISTEQ;
      else if (strncmp(optarg, "lupton", strlen("lupton")) == 0) {
	if (strlen(optarg)>6) {
	  c_lupton = atof (&optarg[6]);
	}
	pscale = LUPTON;
	if (verbose)
	  printf ("Lupton scaling %f\n",c_lupton);
      }
      else if (strncmp(optarg, "pow", strlen("pow")) == 0) {
	if (strlen(optarg)>3) {
	  powpar = atof (&optarg[3]);
	  pscale = POW;
	}
	else
	  printf ("Cannot get value for powpar\n");
      }
      else {
	printf ("Unknown scaling: %s\n", optarg);
	printf ("Accepted values: linear, log, alog, sqrt, square, cbrt, cube\n");
	printf ("                 atanh, atanh, sinh, tanh, pow, lupton\n");
	return (-1);
      }
      break;
    case 'p':
      if (strcmp (optarg, "rainbow") == 0)
	ptype = RAINBOW;
      else if (strcmp (optarg, "rainbow2") == 0)
	ptype = RAINBOW2;
      else if (strcmp (optarg, "rainbow3") == 0)
	ptype = RAINBOW3;
      else if (strcmp (optarg, "grayscale") == 0)
	ptype = GRAYSCALE;
      else if (strcmp (optarg, "red") == 0)
	ptype = RED;
      else if (strcmp (optarg, "green") == 0)
	ptype = GREEN;
      else if (strcmp (optarg, "blue") == 0)
	ptype = BLUE;
      else if (strcmp (optarg, "bb") == 0)
	ptype = BB;
      else if (strcmp (optarg, "heat") == 0)
	ptype = HEAT;
      else if (strcmp (optarg, "bb2") == 0)
	ptype = BB2;
      else if (strcmp (optarg, "IRAF") == 0)
	ptype = IRAF;
      else if (strcmp (optarg, "AIPS") == 0)
	ptype = AIPS;
      else {
	printf ("Unknown palette type: %s\n", optarg);
	return (-1);
      }
      break;
    case 't':
      title = optarg;
      break;
    case 'o':
      output = optarg;
      break;
    case 'a':
      smin = atof (optarg);
      break;
    case 'b':
      smax = atof (optarg);
      break;
    case 's':
      pixelscale_x = atof (optarg);
      pixelscale_y = atof (optarg);
      pixelscale = true;
      break;
    case 'z':
//       pixelscale_y = atof (optarg);
//       pixelscale = true;
      window_size = atoi (optarg);
      break;
    case 'l':
      xlabel = optarg;
      break;
    case 'm':
      ylabel = optarg;
      break;
    case 'h':
      Help ();
      return (EXIT_SUCCESS);
    }
  }

  if (smin!=-1e10 && smin==smax) {
    fprintf (stderr, "WARNING: smin == smax\n  using smin = -smax instead\n");
    smin = -smax;
  }

  char *prompt;
  if (output != NULL) {
    prompt = new char[strlen(output)+7];
    strcpy (prompt, output);
    if (strstr (output, ".gif")!= NULL ) 
      strcat (prompt, "/VGIF");
    else if (strstr(output, ".png")!=NULL) {
      strcat (prompt, "/PNG");
    }
    else {
      strcat (prompt, "/VCPS");
      postscript = true;
    }
  }
  else {
    prompt = new char[6];
    strcpy (prompt, "/XWIN");
  }

  if ( pscale == HISTEQ && postscript == false)
    printf ("Histeq option may not work properly for non-postscript output\n");
  if ( pscale == HISTEQ && postscript == true && postscripthack == false)
    printf ("24bit output is probably required for histeq option\n");

  if (window_size>0 && yref1>0.0 && xref1>0.0) {
    sminx = (int)(xref1 - (float)window_size);
    sminy = (int)(yref1 - (float)window_size);
    smaxx = (int)(xref1 + (float)window_size);
    smaxy = (int)(yref1 + (float)window_size);
    if (verbose) {
      printf ("Window_size: %d\nsminx:%d\nsminy%d\nsmaxx:%d\nsmaxy%d\n",window_size,sminx,sminy,smaxx,smaxy);
    }
  }

  if (verbose)
    printf ("Opening pgplot driver\n");
  if (cpgbeg(0, prompt, 1, 1) != 1)
    exit(EXIT_FAILURE);
  delete[] prompt;

  //cpgldev();
  cpgscf (2);

  if (optind < argc ) {
    bool uselua = false;
    int status = 0; 
    fitsfile *fptr;
    int hdutype;

    if ( strstr(argv[optind], ".lua")) {
#ifndef NO_LUA
      lua.ExecFile (argv[optind]);
      uselua = true;
#else
      printf ("Lua not included during compile\n");
#endif
    }
    else if ( !fits_open_file(&fptr, argv[optind], READONLY, &status) ) {
      if (fits_get_hdu_type(fptr, &hdutype, &status) || hdutype != IMAGE_HDU) { 
	fits_close_file (fptr, &status);
	LoadFile (argv[optind], &data);
      }
    
      else {
	LoadFitsFile (fptr, &data);
	fits_close_file (fptr, &status);
      }
    }
    else 
      LoadFile (argv[optind], &data);

    if (trismooth)
      TriangularSmooth (&data);

    if (zoom > 0.0)
      Resample (zoom, &data, true);

    if (flipx) {
      pixelscale_x = -1.0*pixelscale_x;
      xref = xref + sizex*pixelscale_x;
    }

    if (!uselua) {
      Draw();
#ifndef NO_LUA
      if (luastring!=NULL)
	lua.ExecString (luastring);
      if (luafile != NULL)
	lua.ExecFile (luafile);
      else if (optind == argc-2) {
	if (verbose)
	  printf ("Loading lua script %s\n", argv[optind+1]);
	lua.ExecFile (argv[optind+1]);
      }
#endif
      cpgsch (0.55*labelsize);
      cpgslw (linewidth);
      cpgslw (labelthick);
      
      // wedge
      if (wedge) {
	setvp();
	cpgwnad (0.0+xref, -(sizex+1)*pixelscale_x+xref, 
		 -yref, (sizey+1)*pixelscale_y-yref);
	if (!wedgeright)
	  DrawWedge ("BI", 4.0/0.7, 5.0, flmin, flmax, wedgelabel);
	else
	  DrawWedge ("RI", 0.5, 7.0, flmin, flmax, wedgelabel);
      }
    }
  }

  cpgend();

  if (postscript && postscripthack) {
    if (gimage != NULL && bimage != NULL) {
      float *gdata;
      float *bdata;
      
      float rflmin=flmin;
      float rflmax=flmax;
      
      FileAll (gimage, &gdata);
      FileAll (bimage, &bdata);

      if (trismooth) {
	TriangularSmooth (&gdata);
	TriangularSmooth (&bdata);	
      }
      
      if (zoom > 0.0) {
	Resample (zoom, &gdata, false);
	Resample (zoom, &bdata, false);
      }

      if (bflmin < -1e9) bflmin= rflmin;
      if (bflmax < -1e9) bflmax= rflmax;
      if (gflmin < -1e9) gflmin= rflmin;
      if (gflmax < -1e9) gflmax= rflmax;
      
      std::vector<char> * res =
	PSCharArray3 (sizex, sizey, 
		      data,  rflmin, rflmax, 
		      gdata, gflmin, gflmax, 
		      bdata, bflmin, bflmax, 
		      RGB, RGB_val, nRGB);
      InsertPSCharArray (output, res, NULL);
      delete res;
      delete[] bdata;
      delete[] gdata;
    }
    else {
      std::vector<char> * res =
	PSCharArray (data, sizex, sizey, flmin, flmax, 
		     RGB, RGB_val, nRGB);
      InsertPSCharArray (output, res, NULL);
      delete res;
    }
  }

  if (data != NULL)
    delete[] data;


  if (nRGB != 0) {
    delete[] RGB[0];
    delete[] RGB[1];
    delete[] RGB[2];
    delete[] RGB_val;
  }
  
  return EXIT_SUCCESS;
}

void Draw ()
{
  float tr[6] = {xref, -pixelscale_x, 0.0, -yref, 0.0, pixelscale_y};
  if (verbose) {
    printf ("Draw:\n  xref=%f, yref=%f\n", xref, yref);
    printf ("  sizex=%d, sizey=%d\n", sizex, sizey);

  }

  cpgpage ();
  setvp();
  cpgwnad (-0.5*pixelscale_x+xref, -(sizex+0.5)*pixelscale_x+xref, 
	   -yref+0.5*pixelscale_y, (sizey-0.5)*pixelscale_y-yref);
  cpgsitf (0);

  if (smin != -1e10)
    flmin = smin;
  if (smax != -1e10)
    flmax = smax;

//   flmin = flmin-(flmax-flmin)/(float)ar_size;

  if (verbose)
    printf ("Calling palett\n");
  palett (ptype, pscale, contra, bright);

  if (perc)
    FluxPerc (&flmin, &flmax);

  if (verbose)
    printf ("Calling cpgimag\n");
  cpgimag (data, sizex, sizey, 1, sizex, 1, sizey, flmin, flmax, tr);

  if (contour) {
    // LoadContours (contour);
    // if (contour_levels.size() > 0)
    //  cpgcont (data, sizex, sizey, 1, sizex, 1, sizey, 
    //	       &contour_levels[0], -(int)contour_levels.size(), tr);
#ifndef NO_LUA
    lua.ExecString ("Contour ({});\n");
#else
    printf ("Lua support not included during compilation\n");
    printf ("Unable to draw contours\n");
#endif
  }

  cpgqlw (&linewidth);

  // Annotate plot
  cpgslw (labelthick);
  cpgsch (1.0*labelsize);
  if (title!=NULL) 
    cpgmtxt("t", 1.0, 0.5, 0.5, title);
  cpgsch (0.85*labelsize);

  char com1[32];
  char com2[32];
  char com1l[32];
  char com2l[32];
  com1[0] = 0;
  com2[0] = 0;
  com1l[0] = 0;
  com2l[0] = 0;
  
  if (axis) {
    strcat (&com1[0], "bc");
    strcat (&com2[0], "bcv");
  }
  if (ticks) {
    strcat (&com1[0], "ts");
    strcat (&com2[0], "ts");
  }
  
  if (invert) 
    cpgsci (0);
  cpgsch (ticklength);
  cpgslw (tickthick);
  cpgbox (&com1[0], 0.0, 0, &com2[0], 0.0, 0);
  cpgsch (0.85*labelsize);
  cpgslw (labelthick);
  if (invert) 
    cpgsci (1);
  if (numbers)
    if (noxnumbers)
      cpgbox ("", 0.0, 0, "nv", 0.0, 0);
    else if (noynumbers)
      cpgbox ("n", 0.0, 0, "", 0.0, 0);
    else
      cpgbox ("n", 0.0, 0, "nv", 0.0, 0);

  cpgslw (labelthick);
  if (xlabel!= NULL) 
    cpgmtxt ("b", 3.0, 0.5, 0.5, xlabel);

  if (ylabel != NULL)
    cpgmtxt ("l", 3.0, 0.5, 0.5, ylabel);

  //cpgsch (0.75);

  cpgsch (1.0);
}

void palett (palettetype type, palettescaling scale, 
  float cont, float brg)
{
  int ar_size = 301;

  int it1, it2, ic1, ic2;
  cpgqcol (&ic1, &ic2);
  if (verbose)
    printf ("Max colors: %d\n", ic2);
  cpgscir (16, 316);
  cpgqcir (&it1, &it2);
  if (verbose)
    printf ("Color index: %d - %d\n", it1, it2);
  ar_size = it2 - it1 + 1;


  float v_null[ar_size], v_linear[ar_size], v_linear1[ar_size], v_linear2[ar_size], v_linear3[ar_size];
  float v_cbrt[ar_size], v_p3[ar_size];
  float v_sqrt[ar_size];
  float v_p2[ar_size];
  float v_alog[ar_size];
  float v_log[ar_size];
  float v_pow[ar_size];
  float v_asinh[ar_size];
  float v_sinh[ar_size];
  float v_tanh[ar_size];
  float v_atanh[ar_size];
  float v_lupton[ar_size];
  float v_histeq[ar_size];

  float chr[ar_size];
  float chg[ar_size];
  float chb[ar_size];

  float BB_r[ar_size];
  float BB_g[ar_size];
  float BB_b[ar_size];

  float heat_r[ar_size];
  float heat_g[ar_size];
  float heat_b[ar_size];
  
  float crr[ar_size], crg[ar_size], crb[ar_size];
  float crr2[ar_size], crg2[ar_size], crb2[ar_size];
  float crr3[ar_size], crg3[ar_size], crb3[ar_size];

  float wl[10] = {0.0, 0.5, 0.5, 0.7, 0.7, 0.85, 0.85, 0.95, 0.95, 1.0};
  float wr[10] = {0.0, 1.0, 0.0, 0.0, 0.3,  0.8,  0.3,  1.0,  1.0, 1.0};
  float wg[10] = {0.0, 0.5, 0.4, 1.0, 0.0,  0.0,  0.2,  0.7,  1.0, 1.0};
  float wb[10] = {0.0, 0.0, 0.0, 0.0, 0.4,  1.0,  0.0,  0.0, 0.95, 1.0};

  float al[20] = {0.0, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3, 0.4, 0.4, 0.5,
		  0.5, 0.6, 0.6, 0.7, 0.7, 0.8, 0.8, 0.9, 0.9, 1.0};
  float ar[20] = {0.0, 0.0, 0.3, 0.3, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0,
		  0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
  float ag[20] = {0.0, 0.0, 0.3, 0.3, 0.0, 0.0, 0.0, 0.0, 0.8, 0.8,
		  0.6, 0.6, 1.0, 1.0, 1.0, 1.0, 0.8, 0.8, 0.0, 0.0};
  float ab[20] = {0.0, 0.0, 0.3, 0.3, 0.7, 0.7, 0.7, 0.7, 0.9, 0.9,
		  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

  for (int i=0; i<ar_size; i++) {
    //histeq[i] = 0;
    float val = (float)((1.0/(float)(ar_size-1))*(float)i);
    v_linear[i]  = val;
    v_linear1[i] = val;
    v_linear2[i] = val;
    v_linear3[i] = val;
    v_sqrt[i] = sqrt(val);
    v_asinh[i] = asinhf (val);
    v_sinh[i] = sinhf (val);
    v_tanh[i] = tanhf(10.0*val);
    v_atanh[i] = 0.05*atanhf(val);
    //    v_lupton[i] = asinhf(sinhf(c_lupton)*val)/c_lupton;
    v_lupton[i] = sinhf(c_lupton*val)/sinh(c_lupton);
    v_cbrt[i] = pow (val, 1.0/3.0);
    v_p3[i] = pow (val, 3.0);
    v_p2[i] = pow(val, 2.0);
    v_null[i] = 0.0;
    v_pow[i] = pow(val, powpar);

    v_log[i] = (log10(lgcon+val)-log10(lgcon)) 
      / (log10(1.0+lgcon)-log10(lgcon));
    v_alog[i] = pow(10, (log10(1.0+lgcon)-log10(lgcon))*val+log10(lgcon))-lgcon;

    if (val<=0.39)
      crr2[i] = 0.0;
    else if (val<=0.50)
      crr2[i] = (val-0.39)/(0.5-0.39)*0.3;
    else if (val<=0.6)
      crr2[i] = (val-0.50)/(0.6-0.5)*0.7 + 0.3;
    else
      crr2[i] = 1.0;

    if (val<=0.25)
      crr3[i] = 0.4-1.6*val;
    else if (val<=0.39)
      crr3[i] = 0.0;
    else if (val<=0.50)
      crr3[i] = (val-0.39)/(0.5-0.39)*0.3;
    else if (val<=0.6)
      crr3[i] = (val-0.50)/(0.6-0.5)*0.7 + 0.3;
    else
      crr3[i] = 1.0;
    crr3[0] = 0.0;

    if (val<=0.17)
      crg2[i] = 0.0;
    else if (val <= 0.39)
      crg2[i] = (val-0.17)/(0.39-0.17)*0.8;
    else if (val <= 0.5)
      crg2[i] = (val-0.39)/(0.5-0.39)*0.2+0.8;
    else if (val <= 0.61)
      crg2[i] = 1.0;
    else if (val <= 0.75)
      crg2[i] = 1.0 - (val - 0.61)/(0.75-0.61)*0.4;
    else if (val<=0.85)
      crg2[i] = 0.6 - (val-0.75)/(0.85-0.75)*0.6;
    else
      crg2[i] =  (val-0.85)/(1.0-0.85);

    if (val<=0.17)
      crg3[i] = 0.0;
    else if (val <= 0.39)
      crg3[i] = (val-0.17)/(0.39-0.17)*0.8;
    else if (val <= 0.5)
      crg3[i] = (val-0.39)/(0.5-0.39)*0.2+0.8;
    else if (val <= 0.61)
      crg3[i] = 1.0;
    else if (val <= 0.75)
      crg3[i] = 1.0 - (val - 0.61)/(0.75-0.61)*0.4;
    else if (val<=0.85)
      crg3[i] = 0.6 - (val-0.75)/(0.85-0.75)*0.6;
    else
      crg3[i] =  (val-0.85)/(1.0-0.85);

    if (val<=0.17)
      crb2[i] = val/0.17*0.5 + 0.3;
    else if (val<=0.39)
      crb2[i] = (val-0.17)/(0.39-0.17)*0.2 + 0.8;
    else if (val<=0.5)
      crb2[i] = 1.0 - (val-0.39)/(0.5-0.39)*0.7;
    else if (val<=0.61)
      crb2[i] = 0.3 - (val-0.5)/(0.61-0.5)*0.3;
    else if (val<=0.85)
      crb2[i] = 0;
    else
      crb2[i] = (val-0.85)/(1.0-0.85);

    if (val<=0.17)
      crb3[i] = val/0.17*0.5 + 0.3;
    else if (val<=0.35)
      crb3[i] = (val-0.17)/(0.35-0.17)*0.2 + 0.8;
    else if (val<=0.5)
      crb3[i] = 1.0 - (val-0.35)/(0.5-0.35)*0.7;
    else if (val<=0.61)
      crb3[i] = 0.3 - (val-0.5)/(0.61-0.5)*0.3;
    else if (val<=0.85)
      crb3[i] = 0;
    else
      crb3[i] = (val-0.85)/(1.0-0.85);

    if (val<=0.33)
      crr[i] = 0.0;
    else if (val<=0.50)
      crr[i] = (val-0.33)/(0.5-0.33)*0.3;
    else if (val<=0.67)
      crr[i] = (val-0.50)/(0.67-0.5)*0.7 + 0.3;
    else
      crr[i] = 1.0;

    if (val<=0.17)
      crg[i] = 0.0;
    else if (val <= 0.33)
      crg[i] = (val-0.17)/(0.33-0.17)*0.8;
    else if (val <= 0.5)
      crg[i] = (val-0.33)/(0.5-0.33)*0.2+0.8;
    else if (val <= 0.67)
      crg[i] = 1.0;
    else if (val <= 0.83)
      crg[i] = 1.0 - (val - 0.67)/(0.83-0.67)*0.4;
    else 
      crg[i] = 0.6 - (val-0.83)/(1.0-0.83)*0.6;

    if (val<=0.17)
      crb[i] = val/0.17*0.5 + 0.3;
    else if (val<=0.33)
      crb[i] = (val-0.17)/(0.33-0.17)*0.2 + 0.8;
    else if (val<=0.5)
      crb[i] = 1.0 - (val-0.33)/(0.5-0.33)*0.7;
    else if (val<=0.67)
      crb[i] = 0.3 - (val-0.5)/(0.67-0.5)*0.3;
    else 
      crb[i] = 0;

    if (val<=0.4)
      chr[i] = val*2.5;
    else
      chr[i] = 1.0;      

    if (val<=0.2)
      chg[i] = 0.0;
    else if (val<=0.6)
      chg[i] = val*2.5 - 0.5;
    else
      chg[i] = 1.0;

    if (val<=0.4)
      chb[i] = 0.0;
    else if (val<=0.6)
      chb[i] = val*1.5 - 0.6;
    else 
      chb[i] = (val-0.6)*0.7/0.4 + 0.3;

    if (val<=0.5)
      BB_r[i] = val*2.0;
    else 
      BB_r[i] = 1.0;

    if (val<=0.25)
      BB_g[i] = 0.0;
    else if (val<0.75)
      BB_g[i] = (val-0.25)*2.0;
    else
      BB_g[i] = 1.0;

    if (val<=0.5)
      BB_b[i] = 0.0;
    else
      BB_b[i] = (val-0.5)*2.0;

    if (val<=0.34)
      heat_r[i] = val/0.34;
    else
      heat_r[i] = 1.0;

    heat_g[i] = val;
    
    if (val<0.65)
      heat_b[i] = 0.0;
    else if (val<=0.98)
      heat_b[i] = (val-0.65)/(0.98-0.65);
    else 
      heat_b[i] = 1.0;
  }

  float * fv = NULL;

  switch (scale) {
  case LINEAR:
    fv = v_linear; break;
  case LOG:
    fv = v_alog; break;
  case ALOG:
    fv = v_log; break;
  case SQRT:
    fv = v_p2; break;
  case SQUARE:
    fv = v_sqrt; break;
  case CBRT:
    fv = v_p3; break;
  case CUBE:
    fv = v_cbrt; break;
  case HISTEQ:
    fv = v_histeq; break;
  case POW:
    fv = v_pow; break;
  case ASINH:
    fv = v_asinh; break;
  case SINH:
    fv = v_sinh; break;
  case TANH:
    fv = v_tanh; break;
  case ATANH:
    fv = v_atanh; break;
  case LUPTON:
    fv = v_lupton; break;
  }

  if (scale == HISTEQ) {
    std::vector<float> dump;
    for (int i=0; i<sizex*sizey; i++) {
      float val = (data[i]-flmin)/(flmax-flmin);
      if (val>=0.0 && val<= 1.0) 
	dump.push_back (val);
    }
    sort (dump.begin(), dump.end());

    float tbsize = ((float)(dump.size())) / ((float)ar_size);
    for (int i=1; i<ar_size-1; i++) 
      v_histeq[i] = dump[(int)((float)i*tbsize)];
    v_histeq[0] = 0.0;
    v_histeq[ar_size-1] = 1.0;

 //    int n = 0;
//     int histeq[ar_size];
//     for (unsigned int i=0; i<ar_size; i++) 
//       histeq[i] = 0;
//     for (int i = 0; i<sizex*sizey; i++) {
//       float val = (data[i]-flmin)/(flmax-flmin);
//       if (val>=0.0 && val<= 1.0) {
// 	n++;
// 	histeq[(int)(val*(float)ar_size)]++;
//       }
//     }
//     int sum = 0;
//     int index = 0, index2=0;
//     for (unsigned int i = 0; i<ar_size; i++) {
//       sum = sum + histeq[i];
//       v_histeq[i] = (float)sum/(float)n;
//       //printf ("%d %g\n", i, histeq[i], v_histeq[i]);
//     }
    
//     while (sum!=n) {
//       if (sum<=index2*n/100) {
// 	v_histeq[index2] = (float)sum/(float)n;
// 	index2++;
//       }
//       sum = sum + histeq[index++];
//     }
//     for (unsigned int i = 0; i<ar_size; i++) {
//       sum = sum + histeq[i];
//       v_histeq[i] = (float)sum/(float)n;
//       //printf ("%d %g\n", i, histeq[i], v_histeq[i]);
//     }
  }

  float *val_r, *val_g, *val_b;
  int num = 0;

  switch (type) {
  case GRAYSCALE:
    val_r = v_linear1; 
    val_g = v_linear2; 
    val_b = v_linear3;
    num = ar_size;
    break;
  case HEAT:
    val_r = heat_r;
    val_g = heat_g; 
    val_b = heat_b;
    num = ar_size;
    break;
  case RAINBOW:    
    val_r = crr;
    val_g = crg;
    val_b = crb;
    num = ar_size;
    break;
  case RAINBOW2:    
    val_r = crr2;
    val_g = crg2;
    val_b = crb2;
    num = ar_size;
    break;
  case RAINBOW3:    
    val_r = crr3;
    val_g = crg3;
    val_b = crb3;
    num = ar_size;
    break;
  case RED:
    val_r = v_linear1;
    val_g = NULL;
    val_b = NULL;
    num = ar_size;
    break;
  case GREEN:
    val_r = NULL;
    val_g = v_linear1;
    val_b = NULL;
    num = ar_size;
    break;
  case BLUE:
    val_r = NULL;
    val_g = NULL;
    val_b = v_linear1;
    num = ar_size;
    break;
  case BB:
    val_r = BB_r;
    val_g = BB_g;
    val_b = BB_b;
    num = ar_size;
    break;
  case BB2:
    val_r = chr;
    val_g = chg;
    val_b = chb;
    num = ar_size;
    break;
  case AIPS:
    fv = al;
    val_r = ar;
    val_g = ag;
    val_b = ab;
    num = 20;
    break;
  case IRAF:
    fv = wl;
    val_r = wr;
    val_g = wg;
    val_b = wb;
    num = 9;
    break;
  default:
    printf ("Palettetype not implemented\n");
    abort ();
    break;
  }

  if (reverse) {
    float temp;
    for (int i = 0; i<num/2; i++) {
      if (val_r != NULL) {
	temp = val_r[i];
	val_r[i] = val_r[num-i-1];
	val_r[num-i-1] = temp;
      }
      if (val_g != NULL) {
	temp = val_g[i];
	val_g[i] = val_g[num-i-1];
	val_g[num-i-1] = temp;
      }
      if (val_b != NULL) {
	temp = val_b[i];
	val_b[i] = val_b[num-i-1];
	val_b[num-i-1] = temp;
      }
    }
  }

  if (val_r == NULL)
    val_r = v_null;
  if (val_g == NULL)
    val_g = v_null;
  if (val_b == NULL)
    val_b = v_null;

  if (gamm!=1.0) {
    for (int i=0; i<num; i++) {
      val_r[i] = pow(val_r[i], gamm);
      val_g[i] = pow(val_g[i], gamm);
      val_b[i] = pow(val_b[i], gamm);
    }
  }

  if (whitebcg) {
    val_r[0] = 0.1294;
    val_g[0] = 0.1843;
    val_b[0] = 0.7098;
    val_r[0] = 1;
    val_g[0] = 1;
    val_b[0] = 1;
  }

   if (setbcgcolor) {
     val_r[0] = bcg_r;
     val_g[0] = bcg_g;
     val_b[0] = bcg_b;
   }

   cpgctab (fv, val_r, val_g, val_b, num, cont, brg);

   RGB[0] = new float[num];
   RGB[1] = new float[num];
   RGB[2] = new float[num];
   RGB_val= new float[num];
   for (int i=0; i<num; i++) {
     RGB_val[i] = fv[i];
     (RGB[0])[i] = val_r[i];
     (RGB[1])[i] = val_g[i];
     (RGB[2])[i] = val_b[i];
   }
   nRGB = num;

}


void setvp ()
{
  float d, vpx1, vpx2, vpy1, vpy2;

  cpgsvp (0.0, 1.0, 0.0, 1.0);
  cpgqvp (1, &vpx1, &vpx2, &vpy1, &vpy2);

  d = (vpx2-vpx1<vpy2-vpy1?vpx2-vpx1:vpy2-vpy1)/40.0;

  vpx1 = vpx1 + 5.0*d;
//   if (!wedgeright)
//     vpx2 = vpx2 - 2.0*d;
//   else
    vpx2 = vpx2 - 5.0*d;
  vpy1 = vpy1 + 8.0*d;
  vpy2 = vpy2 - 2.0*d;

  cpgvsiz(vpx1, vpx2, vpy1, vpy2);

}

void LoadFitsFile (fitsfile *fp, float **dd)
{
  int naxis, status = 0;
  long naxes[3],fpixel[3];
  int hdutype, hdunum;

  naxes[0] = 1;
  naxes[1] = 1;
  naxes[2] = 1;

  if (verbose)
    printf ("Loading fits file\n");

  fits_get_img_dim(fp, &naxis, &status);

  fits_get_img_size(fp, naxis, &naxes[0], &status);

  if (naxis == 0) {
    fits_get_num_hdus (fp, &hdunum, &status);
    if (hdunum<2 || status) {
      printf("Error: HDUNUM = %d.\n", hdunum);
      abort ();
    }

    fits_movabs_hdu (fp, 2, &hdutype, &status);
    fits_get_img_dim(fp, &naxis, &status);
    fits_get_img_size(fp, naxis, &naxes[0], &status);
  }

//   if (status || naxis != 2) { 
//     printf("Error: NAXIS = %d.  Only 2D images are supported.\n", naxis);
//     abort ();
//   }

  if (status==1) {
    printf ("ERROR: CFitsio has returned an error in LoadFitsFile.\n"
	    "       Perhaps the file is corrupted?\n");
    abort();
  }

  printf ("  %ld %ld %ld\n", naxes[0], naxes[1], naxes[2]);

  if (sminx == -65536)
    sminx = 1;
  if (sminy == -65536)
    sminy = 1;
  if (smaxx == -65536)
    smaxx = naxes[0];
  if (smaxy == -65536)
    smaxy = naxes[1];

  // Sanity checks
  if (sminx < 1 || sminx>naxes[0])
    sminx = 1;
  if (sminy < 1 || sminy>naxes[1])
    sminy = 1;
  if (smaxx < 1 || smaxx>naxes[0])
    smaxx = naxes[0];
  if (smaxy < 1 || smaxy>naxes[1])
    smaxy = naxes[1];
  if (sminx>=smaxx || sminy>=smaxy ) {
    printf ("Illegal figure x,y limits: using defaults\n");
    sminx = 1;
    sminy = 1;
    smaxx = naxes[0];
    smaxy = naxes[1];
  }

  sizex = smaxx-sminx+1;
  sizey = smaxy-sminy+1;
  //  xref = xref-sminx+1;
  //  yref = yref-sminy+1;

  if (verbose) 
    printf ("  determined: sizex=%d, sizey=%d\n", sizex, sizey);

  *dd = new float[sizex*sizey+10000];
  fpixel[0] = sminx; fpixel[1] = sminy; fpixel[2]= 1;

  flmax = -1e10; flmin=1e10;

  fpixel[0] = 1; fpixel[1] = 1; fpixel[2]= 1;
  float *pix = new float[naxes[0]*naxes[1]];
  fits_read_pix (fp, TFLOAT, fpixel, naxes[0]*naxes[1], 0, pix, 0, &status);
  if (status==1) {
    printf ("ERROR: fits_read_pix has returned an error in LoadFitsFile.\n"
	    "       Perhaps the file is corrupted?\n");
    abort();
  }
  for (int i=sminx; i<=smaxx; i++)
    for (int j=sminy; j<=smaxy; j++) {
      (*dd)[ i-sminx + sizex*(j-sminy)] = pix[i-1+(j-1)*naxes[0]];
    }
  delete[] pix;

  if (smin == -1e10 ||  smax == -1e10) 
    for (int i=0; i<sizex*sizey; i++) {
      if (flmin>(*dd)[i]) flmin=(*dd)[i];
      if (flmax<(*dd)[i]) flmax=(*dd)[i];
    }

  if (verbose)
    printf ("Loading fits file completed\n");

  if (pixelscale) {
    if (xref1 != 0.0 && yref1 != 0.0) {
      // printf ("hahah %d %d\n", sminx, sminy);
      xref = (float)(xref1+1.0-(float)sminx)*pixelscale_x;
      yref = (float)(yref1+1.0-(float)sminy)*pixelscale_y;
    }
    else {
      xref = 0.5*(float)sizex*pixelscale_x;
      yref = 0.5*(float)sizey*pixelscale_y;
    }
  }
  else {
    xref = 1.0;
    yref = 1.0;
    pixelscale_x = -1.0;
  }

  xmin = 1;
  ymin = 1;
  xmax = sizex;
  ymax = sizey;

  if (sigma>0.0 && region>1) {
    if (verbose)
      printf ("  running sigma&region loop\n");
    for (int i=0; i<(sizex-region); i=i+region)
      for (int j=0; j<(sizey-region); j=j+region) {
	for (int it = 0; it<iter; it++) {
	  float sum = 0.0;
	  float sum2 = 0.0;
	  float mean, rms;
	  for (int ii=0; ii<region; ii++)
	    for (int jj=0; jj<region; jj++) {
	      sum = sum+(*dd)[ i+ii + sizex*(j+jj)];
	      sum2= sum2+pow((*dd)[ i+ii + sizex*(j+jj)], 2.0);
	    }
	  mean = sum/pow((float)region, 2.0);
	  rms = sqrt(sum2/pow((float)region, 2.0) -pow(mean, 2.0));
	  for (int ii=0; ii<region; ii++)
	    for (int jj=0; jj<region; jj++) {
	      float temp = fabs(((*dd)[ i+ii + sizex*(j+jj)] - mean)/rms);
	      if ( temp>=sigma)
		(*dd)[ i+ii + sizex*(j+jj)] = mean;
	    }
	}
     }
  }
  
  if (smin == -1e10 && smax == -1e10) 
    printf ("Minimum %g, maximum %g\n", flmin, flmax); 

  for (int i=0; i<sizex; i++)
    for (int j=0; j<sizey; j++) 
      if ( isnan ((*dd)[ i + sizex*j]))
	(*dd)[ i + sizex*j]= flmin-100;

  if (verbose)
    printf ("Loading fits file done\n");
}

void LoadFile (const char *c, float **dd)
{
  std::ifstream is(c);
  std::vector<int> x;
  std::vector<int> y;
  std::vector<float> f;

  if (is == NULL) {
    fprintf (stderr, "Cannot open file %s\n", c);
    exit (-1);
  }

  float  tmp1f, tmp2f, tmp3;

  while ( (is >> tmp1f) && (is >> tmp2f) && (is >> tmp3)  ) {
    x.push_back ((int)(tmp1f+0.1));
    y.push_back ((int)(tmp2f+0.1));
    f.push_back (tmp3);
  }
  is.close ();

  for (unsigned int i = 0; i<x.size(); i++) {
    if (xmin>x[i]) xmin = x[i];
    if (xmax<x[i]) xmax = x[i];
    if (ymin>y[i]) ymin = y[i];
    if (ymax<y[i]) ymax = y[i];
    if (flmin>f[i]) flmin = f[i];
    if (flmax<f[i]) flmax = f[i];
  }

  if (smaxx > -65536 && smaxy > -65536 && sminx > -65536 && sminy > -65536) {
    sizex = smaxx-sminx+1;
    sizey = smaxy-sminy+1;
    xmin = sminx;
    ymin = sminy;
    xmax = smaxx;
    ymax = smaxy;
    if (verbose)
      std::cout << "Limiting 2:" << sizex << " " << sizey<<std::endl;
  }
  else if (smaxx > -65536 && smaxy > -65536) {
    sizex = smaxx;
    sizey = smaxy;
    xmax = smaxx;
    ymax = smaxy;
    xmin = 1;
    ymin = 1;
    if (verbose)
      std::cout << "Limiting 1" << std::endl;
  }
  else {
    sizex = xmax - xmin+1;
    sizey = ymax - ymin+1;
    if (verbose)
      std::cout << "No limiting" << std::endl;
  }

  *dd = new float[ sizex*sizey ];    

  for (int i = 0; i<sizex*sizey; i++)
    (*dd)[i] = -1e30;

  for (unsigned int i = 0; i<x.size(); i++) {
    if (x[i]>=xmin & x[i]<=xmax & y[i]>=ymin & y[i]<=ymax)
      (*dd)[ (x[i]-xmin) + sizex*(y[i]-ymin)] = f[i];
  }

  if (pixelscale) {
    if (xref1 != 0.0 && yref1 != 0.0) {
      xref = (float)(xref1+1.0-(float)xmin)*pixelscale_x;
      yref = (float)(yref1+1.0-(float)ymin)*pixelscale_y;
    }
    else {
      xref = 0.5*(float)sizex*pixelscale_x;
      yref = 0.5*(float)sizex*pixelscale_y;
    }
  }
  else {
    xref = 1.0;
    yref = 1.0;
    pixelscale_x = -1.0;
  }

  if (smin == -1e10 && smax == -1e10) 
    printf ("Minimum %g, maximum %g\n", flmin, flmax); 
//   else {
//     for (int i =0; i<sizex; i++)
//       for (int j =0; j<sizey; j++) 
// 	if ( (*dd)[i+sizex*j]<smin)
// 	  (*dd)[i+sizex*j] =smin;
//   }
}

void Help ()
{
#ifndef NO_LUA
  std::cout << "\nUSAGE: colormap [options] file [luafile]\n\n"
	    << "Please note, that the options are case sensitive\n"
	    << "For a list of commands available in luafile, use command ListCommands();\n\n"
#else
    std::cout << "\nUSAGE: colormap [options] file\n\n"
	    << "Please note, that the options are case sensitive\n"
#endif
	    << "The available command line options are:"
	    << "\n -a <value>  lowest value to display"
	    << "\n -b <value>  highest value to display"
	    << "\n -o <file>   output to file (.gif or .eps/.ps) "

	    << "\n -f <flag>   set flags. Currently supported flags are:"
	    << "\n               flipx:      flip x axis"
	    << "\n               noaxis:     disable displaying x and y axis"
	    << "\n               noaxis:     disable displaying x and y axis"
	    << "\n               nonumbers:  disable displaying the numbers"
	    << "\n                 noxnumbers, noynumbers"
	    << "\n               noticks:    disable drawing ticks"
	    << "\n               nowedge:    disable displaying the color wedge"
	    << "\n               reverse:    reverse colormap"
	    << "\n               rticks:     reverse tickmark color" 
	    << "\n               trismooth:  use 2D triangular kernel to"
	    << "\n                           smooth the data"
	    << "\n               verbose:    verbose output"
	    << "\n               wedgeright: place the wedge vertically on the right"
	    << "\n               8bit:       generate pgplot postscript"
	    << "\n               24bit:      generate 24 bit postscript (DEFAULT)"
	    << "\n -g <code>   set scaling method, possible values are:"
	    << "\n             linear, log, alog, sqrt, square, cbrt, cube "
    	    << "\n             asinh, sinh, atanh, tanh and histeq"
	    << "\n             pow and lupton have optional parameters"
	    << "\n -l <label>  label for the X-axis"
	    << "\n -m <label>  label for the Y-axis"
	    << "\n -p <code>   set color palette type, possible values are:"
	    << "\n             rainbow (default), grayscale, red, green, blue"
	    << "\n             bb, heat, bb2, iraf, aips, rainbow2, rainbow3(?)"
	    << "\n -s <value>  pixelscale for both x and y"
	    << "\n -t <title>  set title of the figure"
	    << "\n -v <var=va> set variables (e.g. \" -v label_thickness=2.0\")"
	    << "\n               currently supported are:"
	    << "\n               contrast"
	    << "\n               brightness"
	    << "\n               contoursigma (smoothing factor)"
	    << "\n               contourstart (low sigma, default 3.0)"
	    << "\n               contoursampling (sampling factor, default 1.0)"
	    << "\n               tick_length"
	    << "\n               tick_thickness"
	    << "\n               label_size"
	    << "\n               label_thickness"
	    << "\n               x1,y1,x2,y1"
	    << "\n               cx,cy"
	    << "\n               sigma,region,iter"
	    << "\n               scale (resampling factor)"
	    << "\n               bcgcolor (e.f.#C0C000)"
	    << "\n -x <value>  reference x pixel"
	    << "\n -y <value>  reference y pixel"
    //	    << "\n -z <value>  pixelscale for y only"
	    << "\n -h          display this help"
	    << std::endl;
}


void Resample (float sf, float **ar, bool update)
{
  image_t * im_in, * im_out;
  double          lineartrans[6];
  float *temp;
  int sfpad;
  int isf;

  int sx, sy;

  int pad = 6;

  sfpad = (int)(sf*pad);
  isf = (int)(sf+0.1);

  if (update) {
    sx = sizex;
    sy = sizey;
  } 
  else {
    sx = oldsizex;
    sy = oldsizey;
  }

  im_in = image_new (sx+2*pad, sy+2*pad);

  for (int i =0; i<sx; i++)
    for (int j =0; j<sy; j++)
      if ( (*ar)[i+j*sx]<-1e29) {
	float ff = 0.0;
	float nn = 0.0;
	if (i-1>=0 && (*ar)[i-1+j*sx]>-1e29) {
	  ff = ff + (*ar)[i-1+j*sx];
	  nn = nn + 1.0;
	}
	if (i+1<sx && (*ar)[i+1+j*sx]>-1e29) {
	  ff = ff + (*ar)[i+1+j*sx];
	  nn = nn + 1.0;
	}

	if (i-1>=0 && j-1>=0 && (*ar)[i-1+(j-1)*sx]>-1e29) {
	  ff = ff + (*ar)[i-1+(j-1)*sx];
	  nn = nn + 1.0;
	}
	if (j-1>=0 && (*ar)[i+(j-1)*sx]>-1e29) {
	  ff = ff + (*ar)[i+(j-1)*sx];
	  nn = nn + 1.0;
	}
	if (j+1<sx && j-1>=0 && (*ar)[i+1+(j-1)*sx]>-1e29) {
	  ff = ff + (*ar)[i+1+(j-1)*sx];
	  nn = nn + 1.0;
	}

	if (i-1>=0 && j+1<sy && (*ar)[i-1+(j+1)*sx]>-1e29) {
	  ff = ff + (*ar)[i-1+(j+1)*sx];
	  nn = nn + 1.0;
	}
	if (j+1<sy && (*ar)[i+(j+1)*sx]>-1e29) {
	  ff = ff + (*ar)[i+(j+1)*sx];
	  nn = nn + 1.0;
	}
	if (j+1<sx && j+1<sy && (*ar)[i+1+(j+1)*sx]>-1e29) {
	  ff = ff + (*ar)[i+1+(j+1)*sx];
	  nn = nn + 1.0;
	}
	if (nn>0.5)
	  im_in->data[i+pad+(j+pad)*(sx+2*pad)] = ff/nn;
	else
	  im_in->data[i+pad+(j+pad)*(sx+2*pad)] = 0.0;
      }
      else
	im_in->data[i+pad+(j+pad)*(sx+2*pad)] = data[i+j*sx];
  
  for (int k=0; k<pad; k++) {
    for (int i=pad;i<sx+pad; i++) {
      im_in->data[i+k*(sx+2*pad)]=im_in->data[i+pad*(sx+2*pad)];
      im_in->data[i+(sy+2*pad-1-k)*(sx+2*pad)]=im_in->data[i+(sy+pad-1)*(sx+2*pad)];
    }
    
    for (int j=pad;j<sy+pad; j++) {
      im_in->data[k+j*(sx+2*pad)] = im_in->data[pad+j*(sx+2*pad)];
      im_in->data[sx+2*pad-1-k+j*(sx+2*pad)] = im_in->data[sx+pad-1+j*(sx+2*pad)];
    }
  }

  lineartrans[0] = sf ;
  lineartrans[1] = 0.00 ;
  lineartrans[2] = 0.00 ;
  lineartrans[3] = 0.00 ;
  lineartrans[4] = sf ;
  lineartrans[5] = 0.00 ;

  if (verbose)
    printf ("Warping\n");
  im_out = image_warp_linear(im_in, lineartrans, NULL);
  if (verbose)
    printf ("Warping done\n");

  temp = new float[(im_out->lx-2*sfpad)*(im_out->ly-2*sfpad)];

  for (int i =sfpad; i<im_out->lx-sfpad; i++)
    for (int j =sfpad; j<im_out->ly-sf*pad; j++)
      temp[(i-sfpad)+(j-sfpad)*(im_out->lx-2*sfpad)] 
	= im_out->data[i+j*im_out->lx];

  for (int i =0; i<sx; i++)
    for (int j =0; j<sy; j++) 
      if (data[i+j*sx]<-1e29) {
	for (int ii = 0; ii<isf; ii++)
	  for (int jj = 0; jj<isf; jj++) {
	    temp[ii+i*isf+(jj+j*isf)*(im_out->lx-2*sfpad)]
	      = -1e30;
	  }
      }

  for (int i =0; i<sx; i++)
    for (int j =0; j<sy; j++) {
      int nei=0;
      for (int ii=-1; ii<2; ii++)
	for (int jj=-1; jj<2; jj++)
	  if ( (i+ii)>=0 && (i+ii)<sx
	       && (j+jj)>=0 && (j+jj)<sy
	       && (jj!=0 && ii!=0)) 
	    if ( (*ar)[i+ii+(j+jj)*sx]>-1e29)
	      nei++;
      if (nei < 2) {
	for (int ii = 0; ii<isf; ii++)
	  for (int jj = 0; jj<isf; jj++) 
	    temp[ii+i*isf+(jj+j*isf)*(im_out->lx-2*sfpad)]
	      = (*ar)[i+j*sx];
      }
    }
  
  if (update) {
    pixelscale_x = pixelscale_x/sf;
    pixelscale_y = pixelscale_y/sf;
    oldsizex = sizex;
    oldsizey = sizey;
    sizex = im_out->lx-2*sfpad;
    sizey = im_out->ly-2*sfpad;
  }

  delete[] (*ar);
  *ar = temp;

  image_del (im_in);
  image_del (im_out);
}


void DrawWedge (const char *side, float disp, float width, 
		float fg, float bg, const char *label)
{
  //bool pgnoto;
  float wxa, wxb, wya, wyb, xa, xb, ya, yb; // Temporary window coord storage.
  float vxa, vxb, vya, vyb; // Viewport coords of wedge.
  float oldch, newch; // Original and anotation character heights.
  float ndcsiz;       // Size of unit character height (NDC units).
  bool horiz = true;  // True if wedge plotted horizontally.
  bool image=true;    // Use PGIMAG (T) or PGGRAY (F).
  int nside=0, bot=1, top=2, lft=3, rgt=4; //  Symbolic version of SIDE.
  int i;
  float wedwid, wdginc, vwidth, vdisp, xch, ych, labwid, fg1, bg1;
  float txtfrc = 0.6; // Set the fraction of WIDTH used for anotation
  float txtsep = 2.2; // Char separation between numbers and LABEL.

  int wdgpix = 300;
  float wdgarr[wdgpix+1]; // Array to draw wedge in.
                        
  float tr[6] = {0.0,1.0,0.0,0.0,0.0,1.0}; // Define the coordinate-mapping function.
  int oldthick;
  int llabelthick;
  int ltickthick;

  cpgsci (1);
  cpgsls (1);

  // Get a numeric version of SIDE.
  if (side[0]=='B' || side[0]=='b') {
    nside = bot; 
    horiz = true;
  }
  else if (side[0]=='T' || side[0]=='t') {
    nside = top;
    horiz = true;
  }
  else if (side[0]=='L' || side[0]=='l') {
    nside = lft;
    horiz = false;
  }
  else if (side[0]=='R' || side[0]=='r') {
    nside = rgt;
    horiz = false;
  }
  else 
    printf ("HAJOO\n");

  ltickthick = (int)std::max(0.8*tickthick, 1.0);
  llabelthick = (int)std::max(0.6*labelthick, 1.0);

  // Determine which routine to use.

  if (strlen(side)<2)
    image = false;
  else if (side[1]=='I' || side[1]=='i')
    image = true;
  else if (side[1]=='G' || side[1]=='g')
    image = false;
  else
    printf ("OHHOH\n");

  cpgbbuf ();

  // Store the current world and viewport coords and the character height.
  cpgqwin (&wxa, &wxb, &wya, &wyb);
  cpgqvp (0, &xa, &xb, &ya, &yb);
  cpgqch (&oldch);
  cpgqlw (&oldthick);


  // Determine the unit character height in NDC coords.

  cpgsch(1.0);
  cpgqcs(0, &xch, &ych); // 0: normalized device coordinates
  if(horiz) 
    ndcsiz = ych;
  else
    ndcsiz = xch;	

  // Convert 'WIDTH' and 'DISP' into viewport units.

  vwidth = width * ndcsiz * oldch;
  vdisp  = disp * ndcsiz * oldch;

  // Determine the number of character heights required under the wedge.

  labwid = txtsep;
  if(label!=NULL) 
    labwid = labwid + 1.0 ;


  // Determine and set the character height required to fit the wedge
  // anotation text within the area allowed for it.
  newch = txtfrc*vwidth / (labwid*ndcsiz);
  cpgsch(newch);

  // Determine the width of the wedge part of the plot minus the anotation.
  // (NDC units).

  wedwid = vwidth * (1.0-txtfrc);

  // Use these to determine viewport coordinates for the wedge + annotation.
  vxa = xa;
  vxb = xb;
  vya = ya;
  vyb = yb;
  if(nside==bot) {
    vyb = ya - vdisp;
    vya = vyb - wedwid;
  }
  else if(nside==top) {
    vya = yb + vdisp;
    vyb = vya + wedwid;
  }
  else if(nside==lft){
    vxb = xa - vdisp;
    vxa = vxb - wedwid;
  }
  else if(nside==rgt) {
    vxa = xb + vdisp;
    vxb = vxa + wedwid;
  }

  // Set the viewport for the wedge.
  cpgsvp(vxa, vxb, vya, vyb);

  // Swap FG/BG if necessary to get axis direction right.

  fg1 = std::max(fg,bg);
  bg1 = std::min(fg,bg);

  // Create a dummy wedge array to be plotted.

  wdginc = (fg1-bg1)/((float)wdgpix);
  for (i=0;i<=wdgpix; i++) 
    wdgarr[i] = bg1 + wdginc*((float)i);


  // Draw the wedge then change the world coordinates for labelling.


  if (horiz) {
    cpgswin(1.0, (float)wdgpix, 0.9, 1.1);
    if (image) 
      cpgimag(wdgarr, wdgpix,1, 1,wdgpix, 1,1, fg,bg, tr);
    else
      cpggray(wdgarr, wdgpix,1, 1,wdgpix, 1,1, fg,bg, tr);
    cpgswin(bg1,fg1,0.0,1.0);
    if (cont_values.size()>0) {
      int style, color;
      cpgqls (&style);
      cpgqci (&color);
      for (unsigned int i=0; i<cont_values.size(); i++) {
	cpgsls (cont_styles[i]);
	cpgsci (cont_colors[i]);

	int n=2;
	float delx[2] = {cont_values[i], cont_values[i]};
	float dely[2] = {0.0, 1.0};
	cpgline (n, delx, dely);
      }
      cpgsls (style);
      cpgsci (color);
    }
  }
  else {
    cpgswin(0.9, 1.1, 1.0, (float)wdgpix);
    if (image) 
      cpgimag(wdgarr, 1,wdgpix, 1,1, 1,wdgpix, fg,bg, tr);
    else
      cpggray(wdgarr, 1,wdgpix, 1,1, 1,wdgpix, fg,bg, tr);
    cpgswin(0.0, 1.0, bg1, fg1);
  }

  // Draw a labelled frame around the wedge.
  if(nside==bot) {
    cpgslw (ltickthick);
    cpgbox("bcst",0.0,0,"bc",0.0,0);
    cpgslw (llabelthick);
    cpgbox("n",0.0,0,"",0.0,0);
    cpgslw (oldthick);
  }
  else if(nside==top) {
    cpgslw (ltickthick);
    cpgbox("bcst",0.0,0,"bc",0.0,0);
    cpgslw (llabelthick);
    cpgbox("m",0.0,0,"",0.0,0);
    cpgslw (oldthick);
  }
  else if(nside==lft) {
    cpgslw (ltickthick);
    cpgbox("bc",0.0,0,"bcst",0.0,0);
    cpgslw (llabelthick);
    cpgbox("",0.0,0,"n",0.0,0);
    cpgslw (oldthick);
  }
  else if(nside==rgt) {
    cpgslw (ltickthick);
    cpgbox("bc",0.0,0,"bcst",0.0,0);
    cpgslw (llabelthick);
    cpgbox("",0.0,0,"m",0.0,0);
    cpgslw (oldthick);
  }

  // Write the units label.

  if(label!=NULL) 
    cpgmtxt(side,txtsep,1.0,1.0,label);

  // Reset the original viewport and world coordinates.

  cpgsvp(xa,xb,ya,yb);
  cpgswin(wxa,wxb,wya,wyb);
  cpgsch(oldch);
  cpgebuf ();
}

bool ParseColor (const char *s, float *r, float *g, float *b)
{
  int dig[8];

  if (strlen(s)!=7 || s[0]!='#') {
    printf ("Not a legal color triplet (e.g. \"#AA77BB\")\n");
    return false;
  }
  for (unsigned int i=1; i<strlen(s); i++) 
    if ( s[i]>0x29 && s[i]<0x40)
      dig[i] = s[i] & 0x0f;
    else if (s[i]>='a' && s[i]<='f') 
      dig[i] = (s[i] & 0x0f) + 9;
    else if (s[i] >='A' && s[i] <= 'F') 
      dig[i] = (s[i] & 0x0f) + 9;
    else {
      printf ("Not a hexadecimal!\n");
      return false;
    }
  
  *r = (16.0*(float)dig[1] + (float)dig[2])/255.0;
  *g = (16.0*(float)dig[3] + (float)dig[4])/255.0;
  *b = (16.0*(float)dig[5] + (float)dig[6])/255.0;
  return (true);  
  
}

void TriangularSmooth (float **data1)
{
  float kernel[9] = {0.0625, 0.125, 0.0625, 0.125, 0.25, 0.125, 
		     0.0625, 0.125, 0.0625};


  float *temp = new float[sizex*sizey];
  float *temp_sum = new float[sizex*sizey];
  for (int i=0; i<sizex*sizey; i++) {
    temp[i] = 0.0;
    temp_sum[i] = 0.0;
  }


  for (int i=0; i<sizex; i++)
    for (int j=0; j<sizey; j++)   
      if (std::isfinite( (*data1)[i+sizex*j]) && (*data1)[i+sizex*j]>-5e9) {
	for (int ii=-1; ii<=1; ii++)
	  for (int jj=-1; jj<=1; jj++) {
	    if ( (i+ii)>=0 && (i+ii)<sizex 
		 && (j+jj)>=0 && (j+jj)<sizey ) {
	      temp_sum[i+ii+sizex*(j+jj)] =
		temp_sum[i+ii+sizex*(j+jj)] 
		+ kernel[ii+1+3*(jj+1)];
	      temp[i+ii+sizex*(j+jj)] =
		temp[i+ii+sizex*(j+jj)] 
		+ kernel[ii+1+3*(jj+1)]*
		(*data1)[i+sizex*j];
	    }
	  }
      }
  
  for (int i=0; i<sizex*sizey; i++)
    if (temp_sum[i]>0.0)
      temp[i] = temp[i] / temp_sum[i];

  delete[] *data1;
  delete[] temp_sum;
  *data1 = temp;


}

void FileAll (const char* name, float **buf) 
{
  int status = 0; 
  fitsfile *fptr;
  int hdutype;

  if ( !fits_open_file(&fptr, name, READONLY, &status) ) {
    if (fits_get_hdu_type(fptr, &hdutype, &status) || hdutype != IMAGE_HDU) { 
      fits_close_file (fptr, &status);
      LoadFile (name, buf);
    }
    
    else {
      LoadFitsFile (fptr, buf);
      fits_close_file (fptr, &status);
    }
  }
  else 
    LoadFile (name, buf);

}

void FluxPerc (float *amin, float*amax)
{
  std::vector<float> arr;
  for (int i=0; i<sizex; i++)
    for (int j=0; j<sizey; j++) {
      if (!isnan(data[i+j*sizex]) && data[i+j*sizex]!=0.0)
	arr.push_back (data[i+j*sizex]);
    }

  sort (arr.begin(), arr.end());
  *amin = arr[minperc*arr.size()];
  *amax = arr[maxperc*arr.size()];
}
